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The quasi-independent curvilinear coordinate approximation (QUICCA) method [K. Nemeth and 
M. Challacombe, J. Chem. Phys. 121, 2877, (2004)] is extended to the optimization of crystal 
structures. We demonstrate that QUICCA is valid under periodic boundary conditions, enabling 
simultaneous relaxation of the lattice and atomic coordinates, as illustrated by tight optimization of 
polyethylene, hexagonal boron-nitride, a (10,0) carbon-nanotube, hexagonal ice, quartz and sulfur 
at the r-point RPBE/STO-3G level of theory. 



I. INTRODUCTION 

Internal coordinates, involving bond stretches, angle 
bends, torsions and out of plane bends, etc. are now rou- 
tinely used in the optimization of molecular structures 
by most standard quantum chemistry programs. Inter- 
nal coordinates are advantageous for geometry optimiza- 
tion as they exhibit a reduced harmonic and anharmonic 
vibrational couplingii2ii4. This effect allows for larger 
steps during optimization, reducing the number of steps 
by 7-10 times for small to medium sized molecules, rela- 
tive to Cartesian conjugate gradient schemes^. 

Recently, Kudin, Scuseria and Schlegel (KSS)^ and 
Andzelm, King-Smith and Fitzgerald (AKF)i reported 
the first crystal structure optimizations using internal 
coordinates. These authors proposed new ways of build- 
ing Wilson's B matrix* under periodic boundary con- 
ditions. The B-matrix (and its higher order gener- 
alizations), defines the transformation between energy 
derivatives in Cartesian coordinates and those in internal 
coordinatesSiS. The scheme presented by KSS allows the 
simultaneous relaxation of atomic positions and lattice 
parameters, while the method developed by AKF allows 
relaxation of atomic positions with a fixed lattice. 

More recently still. Bucko, Hafner and Angyan (BHA)^ 
presented a more "democratic" formulation of Wilson's B 
matrix for crystals, realizing that changes in periodic in- 
ternal coordinates may have non- vanishing contributions 
to the lattice. The BHA definition of the B-matrix will 
be used throughout this paper. 

In the present article, we extend our recently developed 
optimization algorithm, the Quasi Independent Curvilin- 
ear Coordinate Approach (QUICCA)iSi to the condensed 
phase. QUICCA is based on the idea that optimization 
can be carried out in an uncoupled internal coordinate 
representation, with coupling effects accounted for im- 
plicitly through a weighted fit. So far, QUICCA has 
been simple to implement, robust and efficient for iso- 
lated molecules, with a computational cost that scales 
linearly with system size. However, nothing is yet known 
about the applicability of QUICCA for crystals. In crys- 
tals, coupling may be very different from that encoun- 



tered in isolated molecules, and there may be significant 
effects from changes in the lattice. If QUICCA also works 
well for crystals, then it presents a viable offshoot of gra- 
dient only algorithms for the large scale optimization of 
atomistic systems. 

In the following, we first review the construction of 
Wilson's B matrix for crystals (Sec. Ill A|l . and then in 
Sec, mil we go over the QUICCA algorithm and discuss 
its implementation for periodic systems. In Sec. IIVI re- 
sults of test calculations ranging over a broad class of 
crystalline systems are presented. We discuss these re- 
sults in Sec.0 and then go on to present our conclusions 
in Sec. ED 

II. METHODOLOGY 
A. Wilson's B matrix for crystals 

In this section we briefly review the BHA approach to 
construction of the periodic B-matrix, with a somewhat 
simplified notation. A more detailed discussion can be 
found in Ref. and a background gained from Wilson's 
book, Ref. [I. 

The independent geometrical variables of a crystal are 
the fractional coordinates fk, and the lattice vectors hi. 
The absolute Cartesian position of the fc-th atom, fk, is 
related to the corresponding fractional coordinates by 

Tk = h/fe, (1) 

where h = {hi : /i2 : h^} is the matrix of (coloumn wise) 
Cartesian lattice vectors. 

The periodic B-matrix is naturally coloumn blocked, 
with a fractional coordinate block B{r. ^3„^, and a lat- 
tice coordinate block B^^g. The full B-matrix is then 

B = (B-'' : B'') , a, N, X {3Na + 9) matrix with N, the 
number of internal coordinates and Na the number of 
atoms. Elements of B-^ involve total derivatives of inter- 
nal coordinates with respect to fractionals, 

6,, = ^, k^l,Na, (2) 
dfk 
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where <j)i is the i-th internal coordinate and the B-matrix 
elements are atom-blocked as the three- vector hik- Only 
the atoms k that determine internal coordinate i are non- 
zero, making B-'^ extremely sparse. Conventional ele- 
ments of Wilson's B-matrix, S^, are related to these 



total derivatives by the chain rule: 



bik — 



dfk 

Likewise, elements of B'* are 



dfk 



hi = 



; = 1,3 



(3) 
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where is the Z-th component of the fractional coordi- 
nate corresponding to atom m, and the summation goes 
over the set of all atoms that determine i.e. in 
case of torsions and out of plane bendings, there are four 
atoms in this set. From Eq. it is clear that change in 
the lattice has the potential to change each internal co- 
ordinate, through the fractional and depending on sym- 
metry, so that in general we can expect B'* to be dense. 
Also, note that /,„; can be greater then 1 if atom m is 
not in the central cell. 



B. Internal Coordinate Transformations 

With the B-matrix properly defined, transformation 
of Cartesian coordinates, lattice vectors and their corre- 
sponding gradients into internal coordinates remains. As 
in the case of the B-matrix, the Cartesian gradients are 
partitioned into a fractional component. 
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(5) 



and a lattice component, g^, with 9 entries involving the 
total derivatives the emphasis underscores the fact 

dhi 

that some programs produce partial derivatives with re- 
spect to lattice vectors, and this subtly must be taken 
into account (e.g. see Eq. (9) in Ref. 0). 

With this blocking structure, gradients in internal co- 
ordinates are then defined implicitly by the equation 
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(6) 



which may be solved in linear scaling CPU time, through 
Cholesky factorization of the matrix B*B that enters 
the left-handed pseudo-inverse, followed by forward and 
backward substitution, as described in Ref. 11]. Lin- 
ear scaling factorization of B*B may be achieved, as 
the upper left-hand iNa x 3A"a block is hyper-sparse, 
leading also to a hyper-sparse elimination tree and 
sparse Cholesky factors, similar to the case of isolated 
moleculesii. Note however that the factorization of 
BB* , involved in construction of a right-handed pseudo- 
inverse, has the dense 9 x Ni row uppermost and the 



dense iV^ x 9 coloumn foremost, leading to dense Cholesky 
factors. Thus, the left- and right-handed approach are 
not equivalent for internal coordinate transformations of 
large crystal systems, as they are in the case of gas phase 
molecules, when using sparse linear algebra to achieve a 
reduced computational cost. 

Once predictions are made for improved values of in- 
ternal coordinates, based on the internal coordinates and 
their gradients, new Cartesian coordinates are calcu- 
lated via an iterative back-transformation as described in 
Ref. which also scales linearly with the system-size. 
At each step of the iterative back-transformation, frac- 
tional coordinates and lattice vectors are updated and 
the corresponding atomic Cartesian coordinates are com- 
puted. From these, updated values of the internal coordi- 
nates are found, and the next iteration is started. Besides 
these minor details, the back-transformation is the same 
as that used for isolated molecules. 



III. IMPLEMENTATION 
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FIG. 1: Progression of RHF-MIC/STO-3G gradients for a hy- 
drogen bond coordinate in hexagonal ice, starting from the 
crystal structure—. Numbers label QUICCA optimization 
steps, and the star shows the predicted minimum given by 
a weighted line fit (dashed line) of the previous 6 points. 



A. The QUICCA algorithm for crystals 

Details of the QUasi-Independent Curvilinear Coordi- 
nate Approximation (QUICCA) for geometry optimiza- 
tion of isolated molecules have been given in Ref. [lol| . 
Here we provide a brief overview of the method. 

QUICCA is based on the trends shown by internal co- 
ordinate gradients during geometry optimization. See 
Fig. n for example. These trends can be exploited by 
a weighted curve fit for each internal coordinate gradi- 
ent, allowing an independent, one-dimensional extrapo- 
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lation to zero. The predicted minima for each internal 
coordinate collectively determines one optimization step. 
This method works surprisingly well, but only when the 
weights are chosen to account for coupling effects; when 
coupling is strong, the weights should be small and vica 
versa. Also, the fitting process has an important averag- 
ing effect on coupling that contributes to the success of 
QUICCA. 

The only difference in our current implementation of 
QUICCA, relative to that described in Ref. fl^, is that 
merging the connectivities from recent optimization steps 
is no longer carried out. Omitting this connectivity 
merger does not change the results for Baker's test suite 
as given in Ref. 0| , and does not appear to diminish the 
overall effectiveness of QUICCA. 



B. The periodic internal coordinate system 

In setting up a periodic internal coordinate system, it 
is first important to consider the situation where, during 
optimization, atoms wrapping across cell boundaries lead 
to large (unphysical) jumps in bond lengths, angles, etc. 
This situation is avoided here by employing a minimum 
image criterion to generate a set of Cartesian coordinates 
consistent with a fixed reference geometry. 

Also, because internal coordinates span cell bound- 
aries, it is convenient to work with a 3 x 3 x 3 supercell, 
including the central cell surrounded by its 27 nearest 
neighbors. Even though a smaller replica of 8 cells, with 
lattice indices between and 1 contains all necessary lo- 
cal internal coordinates, we prefer to employ the larger 
supercell, to avoid fragmentation of bonds etc at the cell 
boundaries. Then, all internal coordinates are identified 
in the supercell by means of a recognition algorithm, just 
as for isolated molecules. Finally, internal coordinates 
are discarded that do not involve at least one atom in 
the central cell. 

This procedure produces symmetry equivalent inter- 
nal coordinates, among those internal coordinates that 
cross cell boundaries. In the present implementation 
these equivalent coordinates are not filtered out, since 
their presence has no major effect on the optimization; 
the equivalent coordinates result in exactly the same line 
fit and same predicted minima. 

It is worth noting that an appropriate internal coor- 
dinate recognition scheme is extremely important to the 
success of internal coordinate optimization. Here, we are 
using a still experimental algorithm, which we hope to 
describe in a forthcoming paper. 

C. Treatment of constraints 

In the treatment of constraints, we distinguish between 
soft and hard constraints. Soft constraints approach their 
target value as the optimization proceeds, reaching it at 
convergence. Most internal coordinate constraints are of 



the soft type in our implementation. The application of 
soft constraints is particularly useful in situations where 
it is difficult to construct corresponding Cartesian coor- 
dinates that satisfy the constrained values. Hard con- 
straints are set to their required value at the beginning 
and keep their value during the optimization. Cartesian 
and lattice constraints are hard constraints in the current 
implementation. 

Our method of treating hard constraints is similar to 
Baker's projection scheme^^; columns of the B-matrix 
corresponding to hard constraints arc simply zeroed. 
This zeroing reflects the simple fact that a constrained 
Cartesian coordinate or lattice parameter may not vary 
any internal coordinate. Note, that if the lattice parame- 
ters a, b, c, a, P and 7 are constrained, B'* must be trans- 
formed from a lattice- vector representation into a lattice- 
parameter representation by using the generalized inverse 
of the lattice-parameter Jacobian. This transformation 
results in 6 columns corresponding to the 6 lattice param- 
eters. After zeroing the relevant columns, this portion 
of the B'' matrix is transformed back into the original 
Ni X 9 representation. In addition, this approach guar- 
antees that during the iterative back-transformation no 
displacements occur for constrained coordinates or pa- 
rameters. 

In both cases, it is necessary to project out the 
constraint-space component of the Cartesian gradients, 
so that the internal coordinate gradients remain consis- 
tent with the constrained internal coordinates. As recom- 
mended by Pulay^, this is accomplished via projection; 

9' = P9, (7) 

where P is the purification projector that filters out the 
constraint-space. For Cartesian variables, this together 
the aforementioned zeroing of B-^ is sufficient to enforce a 
hard constraint. For constrained internal coordinate vari- 
ables though, this projection is not entirely sufficient, as 
there are further contaminants that can arise in the left- 
handed pseudo-inverse transformation to internal coor- 
dinates. These contaminants can be rigorously removed 
by introducing a further projection in the transformation 
step, as suggested by Bakeri^.. However, toward conver- 
gence these contaminants disappear, and in practice, we 
find good performance without introducing an additional 
purification step. And so, Eq. is the only purification 
used in the current implementation, for both hard and 
soft constraints. 

While purification of the gradients is sufficient for the 
already satisfied hard constraints, the soft internal co- 
ordinate constraints must be imposed at each geome- 
try step. This is accomplished through setting the con- 
strained and optimized internal coordinates, followed by 
iterative back transformation to Cartesian coordinates. 
This procedure finds a closest fit that satisfies the modi- 
fied internal coordinate system, as described in Refs. 
and prllT^ . 
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D. Implementation 

Crystal QUICCA has been implemented in the Mon- 
doSCF suite of linear scaling quantum chemistry codesi^, 
using FORTRAN-95 and sparse (non-atom-blocked) lin- 
ear algebra. Total energies are computed using existing 
fast methods (TRS4i^ ONX", QCTC and HiCui^), and 
the corresponding lattice forces (total derivatives) are 
calculated analytically, with related methods that will 
be described eventually. Full linear scaling algorithms 
have been used throughout. 

These linear scaling algorithms deliver F-point ener- 
gies and forces only. For the Hartree-Fock (HF) model, 
this corresponds to the minimum image criterion (HF- 
MIC)^^. For small unit cells, these F-point effects typ- 
ically lead to different values for symmetry equivalent 
bond and lattice parameters. These effects decay rapidly 
with system size, and are typically less severe for pure 
DFT than for HF-MIC. 

While not always the most efficient option, backtrack- 
ing has been used in all calculations. Backtracking pro- 
ceeds by reducing the steplength by halves, for up to 
three cycles. After that, QUICCA accepts the higher 
energy and carries on. 

In all calculations, the TIGHT numerical thresholding 
schemeii has been used, targeting a relative error of lD-8 
in the total energy and an absolute error of lD-4 in the 
forces. A single convergence criterion is used. That cri- 
terion is that the maximum magnitude of both atomic 
and lattice vector gradients is less than 5D-4 au at con- 
vergence. 

Atomic units are used throughout. 



IV. RESULTS 



A. Test set 



TABLE I: Optimization results for crystal structures at the 
PBE/ST0-3G level of theory and in the F-point approxima- 
tion, using the QUICCA algorithm. 



Molecule 



Number of Optimum energy 
optimization steps (a.u.) 



polyethylene 8 

boron-nitride 5 

(lO-O)carbon-nanotube 7 

ice 15 

quartz 44 

sulfur 89 



-77.56774 
-78.41368 
-1503.75024 
-301.31784 
-1303.03159 
-12595.53586 



We have developed a periodic test suite in the spirit 
of Baker's gas phase set^-. The periodic test set 
includes 6 different systems: Polyethylene, hexago- 
nal boron-nitride, a (10,0)carbon-nanotube, hexagonal 
icei^, quartaSii and sulfusSi. Most of these structures 
were taken either from the Inorganic Crystal Structures 
Database (ICSD)S, or from Cambridge Crystallographic 
Data Centei"2i and the translationally unique positions 
generated with Mercury-*. Details of the geometries used 
are given in Appendix fXl 

Full, simultaneous relaxation of both the lattice and 
atomic positions have been carried out by means of crys- 
tal QUICCA, described above, in the F-point approxi- 
mation at the RPBE/ST0-3G level of theory. 

Table ^ shows the number of optimization steps and 
the optimal energy for each system. While the first four 
test systems converged quickly, quartz and sulfur took 
substantially longer to reach the optimum. In the case of 
quartz, there is a very large (unphysical) deformation, 
wherein four membered rings are formed during opti- 
mization, due perhaps to a combination of a minimal 
basis and F-point effects. The optimized structure of 
quartz is shown in Fig. [5] In the sulfur crystal, 5*8 rings 
interact via a Van der Waals like interaction, which has 
a very flat potential, making this a challenging test case. 




1. Convergence of the energies 

Convergence of the total energy is shown in Figs. I3I5I 
for ice, quartz and sulfur. 



2. Convergence of the gradients 

Figures [6181 show convergence of the maximum Carte- 
sian gradients on atoms and lattice-vector gradients, 
<7max,i, with Optimization step i. 



FIG. 2: The optimal structure of quartz at the F-point 
RPBE/STO-3G level of theory. Note that the picture con- 
tains 8 neighboring unit cells for a better representation of 
intercell bonds. 



B. Urea 

The experimental structure of urea, solved by Swami- 
nanthan et. aB^, has been used as a benchmark for crys- 
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FIG. 3: Convergence of the energy during the optimization of 
hexagonal ice. 
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FIG. 5: Convergence of the energy during the optimization of 
sulfur. 
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FIG. 4: Convergence of the energy during the optimization of 
quartz. 
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FIG. 6: Convergence of the gradients during the optimization 
of hexagonal ice. 



tal optimization by several groups. Kudin, Scuseria and 
Schlegcl^ implemented an early internal coordinate opti- 
mization scheme in the GAUSSIAN programs, and ap- 
plied it to the optimization of RPBE/3-21G urea with 
internal coordinate constraints. At about the same time, 
Civalleri et. aP^ implemented a Cartesian conjugate gra- 
dient scheme in the CRYSTAL program, and carried out 
careful studies examining the effects of k-space sampling 
and integral tolerances on fixed lattice optimization of 
RHF/ST0-3G urea. 

Here, we make direct contact with these works, Refs.j^ 
|25| . However, because the linear scaling methods used 
by MONDOSCF are F-point only, we employ a 2 x 2 x 2 
supercell, so that we may make approximate numerical 
comparison with the k-space methods. This 2x2x2 
supercell involves 16 urea molecules, Ci (no) symmetry, 
128 atoms in total, and more than 850 redundant inter- 
nal coordinates: the number of internal coordinates used 



by QUICCA fluctuates slightly during optimization. For 
comparison, the work of Civalleri et. aS^ makes use of 
P42im symmetry, involving just 8 variables in the fixed 
lattice optimization of urea. In their relaxation of urea, 
Kudin, Scuseria and Schlegel^ employed a 4 molecule cell 
with 5*4 symmetry, with optimization of the lattice and 
atomic centers, but all dihedral angles constrained, in- 
volving 204 redundant internal coordinates. 

Convergence of the RPBE/3-21G MoNDoSCF calcu- 
lations are shown in Fig. in which a full relaxation of 
lattice and atomic centers has been performed, together 
with the energy difference from Ref. 0. The GAUS- 
SIAN values for this calculation, involving constrained 
dihedrals, are -447.6501595 and -447.6632120 for the 
beginning and ending values of the total energy. These 
(and subsequent) values have been normalized to total 
energy per 2 urea molecules. The corresponding MON- 
DOSCF values are -447.648312 and -447.661578. The 
GAUSSIAN energy difference is -0.01305, while the 
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FIG. 7: Convergence of the gradients during the optimization 
of quartz. 
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FIG. 9: Convergence of the energy during full relaxation of 
RPBE/3-21G urea. The bar gives the energy difference com- 
puted by Kudin, Scuseria and Schlegel^. 
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FIG. 8: Convergence of the gradients during the optimization 
of sulfur. 
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MONDOSCF difference is -0.01326. The GAUSSIAN 
optimization converged in 13 steps, while QUICCA took 
24 steps. 

Convergence of the RHF-MIC/ST0-3G MondoSCF 
calculations are shown in Fig. together with the 
energy difference from Ref. (25j . both corresponding 
to relaxation of atomic centers only. The beginning 
and ending CYRSTAL values for this calculation are 
-442.069368 and -442.084595, respectively. For Mon- 
doSCF, they are -442.069473 and -442.084671. The 
energy differences are - .01523 and - .01520 for CRYS- 
TAL and MondoSCF, respectively. The CRYSTAL 
optimization converged in 15 steps, while QUICCA took 
19 steps. 



V. DISCUSSION 

Overall, the behavior of QUICCA for crystal systems 
is similar to that of gas phase systems. For a well be- 
haved system like ice, convergence is rapid and mono- 
tone. For systems undergoing large rearrangements, such 
as the quartz system, QUICCA takes many more steps, 
but still maintains monotone convergence. For both ill- 
conditioned (floppy) gas phase and crystal systems, such 
as sulfur, convergence is slower, with steps that some- 
times raise the energy, even with 3-step backtracking. 

Large amplitude motions can lead to rapid changes in 
curvature of the potential energy surface. In this case, 
the QUICCA algorithm may offer advantages relative to 
strategies based on BFGS-like updates, which are history 
laden. This is because QUICCA employs a weighting 
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scheme that takes large moves into account and that can 
identify recently introduced trends in just a few steps. 

The problems encountered with floppy systems are by 
no means unique to QUICCA, but plague most gradient 
only internal coordinate schemes. Also, as with other 
schemes, we have found the performance of QUICCA to 
be sensitive to the quality of the internal coordinate sys- 
tem. It is our opinion that the difhculties encountered 
with many floppy systems could be overcome with a bet- 
ter choice of internal coordinate system. 

For floppy systems, the ability to resolve small energy 
differences with limited precision (due to linear scaling 
algorithms) can also be problematic. In particular, with 
the TIGHT option, MONDOSCF tries to deliver a relative 
precision of 8 digits in the total energy, and 4 digits of 
absolute precision in forces. For sulfur, achieving atomic 
forces below 5D-4 corresponds to an energy difference 
of lD-5, demanding a relative precision in the total en- 
ergy of lD-10. Exceeding the limits of the TIGHT energy 
threshold can be seen clearly in Fig. |S1 wherein the en- 
ergy has jumps below lD-4. The observant reader will 
also notice that the 21*^* data point is missing from Fig.|51 
This data point was removed because it was lD-4 below 
the converged value of the total energy, -3581 . 2926, con- 
fusing the log- linear plot. These fluctuations are at the 
targeted energy resolution, and are likely due to changes 
in the adaptive HiCu grid for numerical integration of the 
exchange-correlationii. Nevertheless, reliable structural 
information can still be obtained, as absolute precision 
in the forces is retained with increasing system size, al- 
lowing gradient following algorithms such as QUICCA to 
still find a geometry which satisfies the force convergence 
criteria. 

For the urea calculations, very good agreement was 
found between the MondoSCF calculations and the 
CRYSTAL results, in accordance with our previous ex- 
perience for both pure DFT~ and RHF-MIC modelsiSi. 
Slightly less satisfactory agreement was found between 
the MondoSCF and GAUSSIAN calculations, which 
was probably due to the differences in constraint. In both 
cases, the QUICCA calculations took more steps; 4 more 
than CRYSTAL, and 11 more than GAUSSIAN. It 
should be pointed out though, that the MondoSCF cal- 
culations involved a substantially more complicated po- 
tential energy surface: Firstly, the F-point surface lacks 
the symmetry provided by k-space sampling. Secondly, 
the 2x2x2 calculation has many more degrees of free- 
dom, and in particular, lower frequency modes due to the 
larger cell. 



VI. CONCLUSIONS 

We have implemented the Bucko, Hafner and Angyaii^ 
definition of periodic internal coordinates in conjunction 
with the QUICCA algorithm, and demonstrated efficient, 
full relaxation of systems with one, two and three di- 
mensional periodicity. In general, we have found that 
QUICCA performs with an efficiency comparable to that 
of similarly posed gas phase problems, and speculate that 
further enhancement may be achieved through an im- 
proved choice of internal coordinates. 

We have argued that linear scaling internal coordi- 
nate transformations for crystal systems can be achieved 
with a left-handed pseudo-inverse, as the dense rows and 
columns of the periodic B*B matrix determine just the 
last few pivots of the corresponding Cholesky factor. 

We have carried out supercell calculations using a full 
compliment of linear scaling algorithms, including sparse 
linear algebra, fast force and matrix builds and found 
good agreement with k-space methods, involving a mod- 
est number of optimization cycles. Thus, in addition to 
further demonstrating the stability of our linear scaling 
algorithms, we have established QUICCA as a reliable 
tool for large scale optimization problems in the con- 
densed phase. 

In conclusion, QUICCA is a new gradient only ap- 
proach to internal coordinate optimization that is robust 
and generally applicable, both to gas-phase molecules 
and systems of one, two and three dimensional periodic- 
ity. It allows for flexible optimization protocols, involv- 
ing simultaneous relaxation of lattice and atom centers, 
constrained lattice with relaxation of atom centers, con- 
strained atom centers with optimization of the lattice, 
admixtures of the above with constrained internal coor- 
dinates, etc. QUICCA is conceptually simple and easy 
to implement. Perhaps most importantly though, it is a 
new approach to gradient only internal coordinate opti- 
mization, offering a number of opportunities for further 
development. 
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APPENDIX A: TEST SET COORDINATES 

Here, input geometries for the crystal optimization test 
suite are detailed. These geometries are available as sup- 
plementary data, and are also available from the authors 
upon request. 

a. Polyethylene The 1-D periodic structure of 
polyethylene is given in Table ITll 

b. Hexagonal boron-nitride The 2-D periodic coor- 
dinates for hexagonal boron-nitride are given in Table lTTTI 

c. (10,0)carbon-nanotube The geometry of the 1-D 
periodic (10,0)carbon-nanotube has all bond-lengths par- 
allel to the nanotubc axis initially at 1.480A, while those 
running perpendicular to the axis are 1.402 A long. The 



lattice length is a = 4.44A, with the elementary cell con- 

TABLE IL Atomic coordinates in A for the elementary unit 
cell of polyethylene, with corresponding lattice parameters 
a = 2.0A, fe = c = 0.00, a = /3 = 7 = 90.0°. 



c 


0, 


,500 





.500 


0, 


,000 


H 


0, 


,500 


1 


.300 


0, 


,800 


H 


0, 


,500 


1 


.300 


-0, 


,800 


C 


1, 


,500 


-0 


.500 


0, 


,000 


H 


1, 


,500 


-1 


.300 


0, 


,800 


H 


1, 


,500 


-1 


.300 


-0, 


,800 



TABLE III: Fractional coordinates for the elementary unit 
cell of hexagonal boron-nitride with corresponding lattice pa- 
rameters a = 6 = 2.420A,c = 0.00, a = /3 = 90.0°, and 
7 = 120.0°. 

B 0.33333333333 0.1666666666 0.00 
N 0.66666666666 0.8400000000 0.00 



taining 40 atoms. While this data entirely determines the 
structure of the symmetric (10,0)carbon-nanotube, it is 
also available as supplementary data. 

d. Ice Hexagonal ice is the most important natural 
occurrence of ice. Its structure has been taken from the 
literaturei^. Since the literature provides two equilibrium 
position for each hydrogen atom, due to the tunneling of 
hydrogens in ice, our starting structure is taken as the 
average of these two positions for each hydrogen atom, 
and is given in Table HVl 

TABLE IV: Atomic coordinates in A for the elementary unit 
i/jg^j-g^^^/! agonal ice, with corresponding lattice parameters 
a = b = 4.51lA,c = 7.3461, a = 90.0°, /3 = 90.0° and 7 = 
120.0° 
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,0003594 


3, 


,673 


H 


1, 


.1280 


1, 
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1, 


,1280 


1, 


,9540000 
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,346 
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-1, 


,1280 


1, 


,9540000 
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,346 
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,0000 
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,803 
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,6040000 
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,183 
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,2555 
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,0003594 
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,346 



e. Quartz The initial structure of quartz was taken 
from Ref. [23|, and has 9 atoms in the unit cell. 

/. Sulfur The structure of sulfur was taken from 
Ref. [23|. Containing 32 atoms, this structure is avail- 
able as supplementary data. 
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TABLE V: Atomic coordinates in A for the elementary unit 

cell of quartz, with corresponding lattice parameters a = b = 
4.9019A, c = 5.3988A, a = /3 = 90.0° and 7 = 120.0° . 
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.091 


3, 


.094 


2, 


.427 
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